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FOREWORD 

The work described in this report was performed at the Grumman Aerospace 
Corporation, Bethpage, New York, and administered by the Vibration Section of the 
Structures and Dynamics Division, NASA Langley Research Center, Hampton, 

Vi rginia. 

The work performed under NASA Contract NASl-10635-21 with supplementary 
funding provided by the Space Division, Rockwell International (POM3WXMZ-483002) 
included development of a fundamental finite element hydroelastic formulation appli- 
cable to NASTRAN, implementation of the theoretical developments into NASTRAN, 
and verification and demonstration of the new technique on various problems including 
the 1/8-scale space shuttle external tank model. 
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ABSTRACT 

A finite element hydroelastic analysis formulation is developed on the basis of 
Toupin's complementary variational principle of classical mechanics. Emphasis is 
placed on the special case of an incompressible fluid model which is applicable to 
propellant tank hydroelastic analysis. A concise fluid inertia representation results 
from the assumption of incompressibility and the hydroelastic equations reduce to a 
simplified form associated with non-fluid filled structures. The efficiency of the 
incompressible hydroelastic formulation is enhanced for both fluid and structure by 
introduction of harmonic reduction as an alternative to Guyan reduction. The 
theoretical developments are implemented in NASTRAN and the modified NASTRAN 
hydroelastic analysis technique is verified and demonstrated as an efficient and 
accurate approach with a series of illustrative problems including the 1/8-scale 
space shuttle external tank. 
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1 - INTRODUCTION 

The increasing complexity of launch vehicle configurations, particularly in the 
case of the space shuttle, recently has stimulated considerable interest in the dynamic 
behavior of liquid filled tanks. The task of Pogo prediction and suppression, for ex- 
ample, requires very complete and accurate mathematical models for the calculation 
of propellant tank hydroelastic modes in the Pogo-susceptible frequency range (2 - 50 
Hz for the space shuttle) . 

A variety of automatic fluid modeling techniques has been under development 
ranging from finite element and finite difference techniques to approximate analytical 
approaches taking advantage of the properties of the fluid velocity potential and the 
consequences of Green's theorem (Refs. 1-5). The available hydroelastic analysis 
methods although in most cases theoretically rigorous, contain serious deficiencies in 
computational economy and/or numerical accuracy. For example, the NASTRAN hy- 
droelastic analysis technique as formulated in the level 15 series is deficient in 
computational economy primarily because of an unsymmetrical eigenvalue problem re- 
sulting from the use of mixed pressure and displacement generalized coordinates. 

In the NASTRAN hydroelastic formulation, the fluid coefficient matrices are in- 
terpreted according to a structural analogy. The fluid pseudo-mass and pseudo- 
stiffness matrices of that formulation are recognized herein as flexibility and inverse 
mass matrices, respectively, on the basis of the complementary principle in mech- 
anics known as Toupin's principle (Ref. 6). This revised interpretation is central to the 
formulation of the hydroelastic problem presented here. 

This report consists of a theoretical development, a description of NASTRAN 
program modifications, a program with which some familiarity is assumed, and a se- 
ries of illustrative hydroelastic problems demonstrating the accuracy and efficiency of 
the present formulation. The theoretical sections include a derivation of the NASTRAN 
fluid matrix equations on the basis of Toupin's principle, a symmetric formulation for 
compressible hydroelasticity and a symmetric kinematic formulation for incompressi- 
ble hydroelasticity. The incompressible formulation, particularly applicable in the 
study of propellant tank dynamics, provides a description of fluid inertia in terms of 
bounding surface displacements alone. This description represents a drastic reduction 
in system variables. In addition, harmonic reduction is introduced as an efficient 
alternative to Guyan reduction for geometrically axisymmetric structures to further 
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reduce the number of system variables. Detailed theoretical discussions in the appen- 
dices include a derivation of Toupin's principle and the proposed utilization of polyhe- 
dral heat conduction elements in NASTEAN as fluid elements (according to a heat 
conduction-incompressible flow analogy) to accommodate the analysis of asymmetric 
fluid geometries such as a tilted free surface. 

Volume II (Ref. 7) consists of detailed information pertinent to the 
NASTEAN program. Included is a description of NASTEAN program modifications, 
based on the above theoretical developments, consisting of DMAP modifications re- 
quired in the calculation of fluid matrix data and in the calculation normal modes for 
fluid-filled structures with and without the effects of static pressurization. Special 
input bulk data considerations are discussed as an aid to the NASTEAN user. Bulk 
data listings for the illustrative problems are presented in the appendices and serve 
as supplementary user information. 

A series of illustrative hydroelastic problems are presented in the final sec- 
tions of this report. They have been chosen to verify the reformulated NASTEAN 
hydroelastic analysis and to demonstrate its economy. Exact analytical and available 
test results were used as verification data for the NASTEAN analysis. The relatively 
complex 1/8 -scale space shuttle external tank model is included with the illustrative 
examples in spite of a lack of totally satisfactory correlation with experimental data. 
Correlations with exact analytical results and experimental results for all other illus- 
trative examples, however, are excellent and it is concluded that the formulation of 
this report is an accurate and efficient operational approach. 


1-2 
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2 - THEORETICAL DEVELOPMENT 


The class of problems considered in the NASTRAN hydroelastic analysis tech- 
nique consists of the interaction of irrotational, inviscid, compressible fluids with 
flexible structures for which both fluid and structural motions are assumed small 
compared to overall dimensions. The approach used to describe the dynamics of the 
fluid is a finite element technique. "Mass" and "stiffness" matrices are formed on 
the basis of a constructed energy principle with pressure taking the role of general- 
ized displacement and bounding surface displacement taking the role of the forcing 
function. The dynamics of the structure is described in the usual way with displace- 
ment taken as the dynamic variable and applied pressure taken as the forcing function. 
The assembled set of hydroelastic dynamic equations consists of coupled fluid pressure 
and structural displacement matrix relationships containing unsymmetric coupling 
terms as a result of the mixed set of variables. The unsymmetric form of the 
NASTRAN hydroelastic equations leads to considerable analytical and numerical diffi- 
culty. 

2.1 DERIVATION OF THE NASTRAN FINITE ELEMENT FLUID REPRESENTATION 

The NASTRAN fluid equations are derivable on the basis of a complementary 
variational principle introduced by Toupin in 1952, Ref. 6. The physical interpreta- 
tion of the fluid matrix relationships on the basis of this principle provides the insight 
required to resolve the difficulties present in the NASTRAN formulation. (A detailed 
derivation and discussion of Toupin's principle and its consequences is presented in 
Appendix A) . 

The equation of motion of a fluid particle is 

( 2 . 1 - 1 ) 

^f 

The constitutive relationship for an inviscid, compressible fluid is 

P = -B7*U (2.1-2) 

where V^.U represents the dilitational strain. In order to obtain a fluid velocity ex- 
pression, with P£ taken as a mean fluid density, the equation of motion Eq. 2.1-1 is 
integrated resulting in 

U=-^7P (2.1-3a) 

^f 
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where P is the pressure impulse 


or 


-f 


A 

P = P 


P dt 


(2.1-3b) 


(2.1-3C) 


The usual expressions for kinetic and strain energy may now be expressed in terms of 
impulsive pressure as 


T = T^=|/‘ Pj(U-U)dV=iJ i(PP.pP)dV 

V ^f . 

^ (P)^dV 


(2.1 -4a) 
(2.1 -4b) 


V 'V 

The motion dependent and impulse dependent energy expressions are generally not 
equivalent; they are equivalent, however, for linear systems. The complementary 
virtual work performed by boundary surface displacements, U*, is 

A 4- 

d 

'S 


(jw = r<jp(u* - S ) ds 

C Jg 

The complementary form of Hamilton's principle due to Toupin is 


(2.1-5) 


i L dt + iWcdt = 0 
J t. c J j 

h 


with 


L = T - U 
c c c 


(2.1-6a) 


(2.1-6b) 


dS dt= 0(2. 1-7) 


The expression of the principle in the present application is 

(U*-n)dP 

Upon utilization of Green's theorem, integration by parts, and rearrangement of terms 

A 

the final expression (taking ^ P = 0 @ t = tQ, ti) is 


fh 

rr/i A 

1 r 

/ 1 A A , 

A 1 

^‘0 

1 

1 

p ^P 1 (J P dV + J 
^f ' s 

+ u;) 

3P ds J 


dt = 0 (2.1-8) 


2-2 



Volume I 


aad by setting the integrands to zero the well known field equation and natural boundary 
conditions for an inviscid, incompressible fluid result which are 


±V-^V^P = 0mV (2.1-9a) 

J3 

U* = — ^ VP • C or P prescribed on S (2. l-9b) 
n Pj 


The usefullness of the complementary principle lies in approximate analysis rather 
than in the derivation of field equations. Consider an approximation of a fluid pres- 
sure (impulse) state in terms of a finite set of variables. The fluid complementary 
kinetic and strain energies are the quadratic functions 


1 A A 

I ] 


with the symmetric inertance matrix defined as 

2 , 


L.. = 
ij 


d T 

c 

A A 

dP. dP. 
1 J 


and the symmetric flexibility matrix defined as 




C.. 

ij 




(2. 1-lOa) 
(2. 1-lOb) 


(2. 1-lOc) 


(2. 1-lOd) 


The elements of the inertance matrix are proportional to — , and the elements of the 
flexibility matrix are proportional to • The complementary virtual work is ex- 
pressed as 


<5W 


I 

i 


A 


/ ^ u* 

-'COT) 


SdP 


dS(jP. 


(2. 1-lla) 


For the special case in which the surface displacements are physically discretized the 
complementary virtual work may be expressed as 


k i 


(2. 1-llb) 
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with the generalized area matrix defined as 


A., = 
ik 


I(4](f 




(2. 1-llc) 


Substitution of Eq. 2. 1-10 and Eq. 2. 1-11 into Eq. 2. 1-6 with the appropriate 
integrations by parts results in the complementary Euler- Lagrange equations 


By taking the time derivative of this expression noting Eq. 2. 1-3C, the Euler- Lagrange 
equations become 

This is the form of the fluid dynamic finite element equations for individual ele- 
ments and stacked systems of elements in NASTRAN. In the case of a stacked system 
of elements the matrix Aik represents only bounding surface generalized areas and uf 
represents discrete surface displacements. The pressures Pj comprise the set of 
boundary surface and internal pressures; therefore the matrix A^k is rectangular. The 
physical interpretations of the matrix quantities, however differ from the interpreta- 
tions in the NASTRAN theoretical manual. The Cy matrix is a flexibility matrix and 
the L|j matrix is an inverse mass matrix (see Appendix A). This point is realized with- 
out interpretation of Toupin's principle by examination of two special illustrative 
cases. Consider first static deformation (T^.-^0) in which Eq. 2. 1-I2b takes the form 




(2.1-13a) 

which twice integrated is 

r C P = - £ A UJ (2. l-13b) 

J J J ^ 
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The matrix C^. is recognized as a flexibility matrix in terms of pressure. Fur- 
ther clarification is realized by noting that 




(2. 1-14) 


where is a stiffness matrix relating surface F^ forces and displacements Uj^> The 
second special case consists of an incompressible fluid (Cjj-"0) in which Eq. 2. l-12b 
takes the form 


j ^ k 

The inertance matrix, must contain one singularity. This singularity, which will 
be discussed fully in Section 2.3, is due to the fact that an incompressible fluid under 
uniform pressure does not deform. It is however apparent in Eq. 2, 1-15 that the in- 
ertance matrix, as in the case of the flexibility matrix is an inverse "mass-type" 
matrix. 

2. 2 A SYMMETRIC FORMULATION FOR COMPRESSIBLE HYDROELASTICITY 

The formulation presented in the NASTRAN theoretical manual utilizes the com- 
plementary Euler-Lagrange equations for a fluid 

(Lj) {p}- + (Cj) |p} = - (a'^){u} (2. 2-la) 

and a standard set of structural dynamic equations 

(Mg)|u} + (Kg){u^ = (A) |p| (2. 2-lb) 

The above are combined to form the unsymmetric set of hydroelastic equations 



( 2 . 2 - 2 ) 
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Considerable numerical and analytical difficulty has been encountered by NASTRAN 
users due to the unsymmetric form of these equations. 

An alternate symmetric formulation is derivable by the complementary principle 

or by manipulation of the structural dynamic equations. Taking the latter approach, 

the internal structural generalized forces, F , are related to the structural displace- 

s 

ments, U, according to 

(Kg)|u}={F^| (2.2-3a) 

Suppose that K represents a supported stiffness matrix; the transformation to internal 
s 

forces is therefore defined as 

(2.2-3b) 

Substitution of the above into Eq. 2.2-lb results in 

and premultiplication by the inverse of the structural mass matrix yields 

^ '{“I (2.2-5a) 

with 


(Cs) = (K'^ (2.2-5b) 

representing the structural flexibility matrix and 

(Lg) = (M“^) (2.2-5C) 

representing the structural inertance matrix. Utilizing Eq. 2. 2-3b and Eq. 2. 2-5a the 
expression for acceleration to be substituted into the fluid Eq. 2.2-la is 


2-6 
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and the fluid dynamic equation Is rewritten as 

(Lj){p| + (Cj)|p| = - (a'^L^A){p} + (A’^)(Lg)|F J (2.2-7) 

The set of hydroelastic equations in terms of force type variables consisting of 
Eq. 2. 2-5a and Eq. 2. 2-7 now takes the symmetric form 



The formulation presented here provides a basis for modification of the 
NASTRAN formulation for inviscid, compressible fluid hydroelastic problems. The 
class of problems of interest in the current work, however, is limited to inviscid, 
incompressible fluids interacting with flexible structures and an alternate simplified 
kinematic formulation is derivable for this case. 

2.3 A SYMMETRIC KINEMATIC FORMULATION FOR INCOMPRESSIBLE HYDRO- 
ELASTICITY 

The complementary Euler- Lagrange matrix equation set for the special case of 
an incompressible fluid (C£)-^0) in a conveniently partitioned form is 



The pressure partitions P,, P and P correspond to free surface, structural 
interface surface and internal fluid pressure sets, respectively, and the displacement 
partitions and correspond to the free surface and structural interface surface 
displacement sets, respectively. The structural dynamic equation set with applied 
fluid pressure loading is in partitioned form 


/.'’.l.'.A&.l 

V'>!“s/Ki 


+ 




(2.3-2) 
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where M is the structural mass matrix, K is the structural stiffness matrix and 
^ s s 

Kff, Kjg, are additional stiffness matrix contributions resulting from the gravi- 
tational potential and possibly ullage pressure fluctuation. 

It is obvious from Eq. 2. 3-1 that the internal pressures are related to the sur- 
face pressures as 


( S ) 


(2. 3-3) 


and the reduced fluid dynamic equation set in terms of surface quantities only is 







(L..) ^(L.pL. , 
11 ' if IS) 


(2. 3-4a) 


(2. 3-4b) 


The reduced inertance matrix is singular since an incompressible fluid under 
uniform pressure does not deform. This singularity must occur ideally in the indivi- 
dual finite element inertance matrices as it is a necessary condition for compatibility. 
This is analogous to the condition imposed on structural finite elements that requires 
no internal stresses under rigid body motion (i.e. singular element stiffness matrices). 
There are two approaches to eliminating the uniform pressure singularity; the first is 
a general approach and the second is a special case in which part of the fluid surface 
is at zero pressure. 
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In the special case of zero free surface pressure (gravitational and ullage stiff- 
ness negligible) the fluid surface dynamic equation set Eq. 2. 3-4a is under the single 
point constraints 

(2.3-5) 

and the singular pressure state, uniform pressure is removed from the total surface 
pressure set. The partitions of Eq. 2. 3-4a now reduce to 



= - (A5){Uf} - 

(2. 3-6a) 

<^ss>{Ps) 


(2. 3-6b) 


If the free surface displacements are normal to the free surface then the generalized 
area coupling partitions are null ((A^j) = (0), (A^^) = (0)). The matrix partition, (L^g), 
is not singular due the constraint Eq. 2. 3-5 and may now be inverted in Eq. 2. 3-6b 
to form the structural interface pressure recovery equation set 

Substitution of this relationship into Eq. 2. 3-6a yields 

and since consists of outward normal displacements only the free surface area 

matrix is non-singular; thus the free surface acceleration (or displacement) recovery 
relationship is 


(2. 3-7) 


(2. 3-8a) 



(2. 3-8b) 
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The structural dynamic equation set Eq. 2, 3-2 under the constraints Eq. 2. 3-5 
is 

<“ss»K}*<"ss>K}=<v>K} 

Substitution of the pressure recovery relationship Eq. 2. 3-7 into the above set re- 
sults in 

' {»} (2'3-lOa) 

with the symmetric fluid mass matrix formed as 

<“f> = <2- 3-lOb) 

Thus in the special case of an Incompressible, inviscid fluid with zero pressure 
fluctuation on a free surface, the hydroelastic dynamics problem reduces to a standard 
problem in structural dynamics. 

Consider now the general case of incompressible fluid/structure interaction in 
which the entire bounding surface of the fluid interacts with flexible structure and a 
gravitational potential or ullage volume. As in the above special case, all internal 
pressures are dependent on surface pressures Eq. 2. 3-3 and the resultant reduced 
fluid inertance matrix is singular. In the case of a fluid represented by discrete sur- 
face pressures the normalized uniform pressure state is 



Under such loading the surface normal accelerations must be null and the necessary 
property of the reduced inertance matrix is 


£ L.: = 0 


(2. 3-12) 
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for all i. The fluid representations in NASTRAN are currently limited to axisymme- 
tric geometries in which fluid pressures are described in terms of circumferential 
harmonics, e. g. , 

N N 

P(r.,0.,z.)= ^ P (r ,z ) cosK^. + L P* (r , z AsinK^^ (2.3-13) 

' ' ' K=0,1,.. K ' K=l,2,.. 

The fluid finite elements are rectangular and triangular solids of revolution and the 
generalized pressure degrees of freedom are the amplitudes Pj^ and P^. In this case 
the zeroth harmonic (n = O) generalized pressure subset contains the uniform pressure 
state and the higher harmonic (n> O) subsets by definition do not contain the uniform 
pressure state. Thus the discrete pressure considerations described by Eq. 2. 3-11 and 
Eq, 2. 3-12 hold for the zeroth harmonic. The singularity in the reduced inertance 
matrix is removable by artificially setting one pressure to zero (in the case of gener- 
alized harmonic pressures this pressure must be in the zeroth harmonic subset). For 
convenience let the pressure partitions |Pf| |^g| represent the single artificially 
nulled pressure and the remaining pressures, respectively; in addition let |Pj| rep- 
resent a pressure associated with the fluid free surface. The artificially constrained 
pressure set is denoted as 

(2. 3-14) 



in which the primed set represents pressure deviations from a uniform pressure state 
of value, Pj, i.e. , 


P = P - Pf (2.3-15) 

s. Si f 

T T 

Assuming that the area coupling partitions are null ((A^j) = 0, (Aj^) =(0)) the pressure 
deviations are related to the surface displacements as 


hi 

(2. 3-16a) 

)|p' i = - (a"^ )iu 1 

ss'\ s( ' ss'\ sj 

(2. 3-16b) 


2-11 



Volume I 


Solution for the pressure deviations in Eq. 2. 3- 16b and substitution of the result into 
Eq. 2. 3-16a yields the pressure deviation recovery relationship 

= (2-3-17a) 

and the displacement recovery relationship 





(2.3-17b) 


The above result is equivalent to imposing a kinematic constraint on outward normal 
surface flow (incompressibility) 


U 

ss s 



IS ss 

0 

/ 

/ 


U 

ss s 


The companion pressure constraint expression is 


(2. 3-18) 



(2. 3-19) 


Substitution of Eq. 2.3-18 and Eq. 2. 3-19 into Eq. 2. 3-4a yields the result already 
obtained in Eq. 2.3-16b 



where the singularity condition is expressed as 


(L„)-(Lf;KL;^)-l(L;f)=(0) 


(2. 3-20a) 


(2. 3-20b) 
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The constraints presented in Eq. 2. 3-18 and Eq. 2. 3-19 are now applied to the struc- 
tural dynamic equation set Eq. 2. 3-2 as 



(2. 3-21a) 


where 




(2. 3-21b) 




(2. 3-21c) 


The coefficient matrix on the right side of Eq. 2. 3-21a reduces to 

^ (r A„ + I = (0 I A^^) 


( 2 . 3 - 22 ) 


Substitution of the pressure deviation recovery relationship Eq. 2. 3-17a into Eq. 

2. 3-21a noting Eq. 2. 3-22 results in the symmetric kinematic equation set 

(Mg + M^) |Ug| + {Kp = |o} (2. 3-23a) 

where the fluid mass matrix is as in the special case, Eq. 2. 3-lOb 


(Mf) = A 


ss 



T 

ss 


(2. 3-23b) 


and the hydroelastic stiffness matrix is 

r ^ r \^ + k ^, ( 2 . 3-230 
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Complete displacement recovery is obtained through Eq. 2. 3-17b and pressure devia- 
tion recovery is obtained through Eq. 2. 3~17a. Complete pressure recovery is 
achieved by combining the upper partition equation set in Eq. 2. 3-2 with Eq. 2. 3-17 
and Eq. 2. 3-19; thus the pressure recovery equation set is 


I Pf I Kjj r+ |u^ j. 



(2. 3-24a) 


(2. 3-24b) 


The general symmetric kinematic formulation developed above is useful in hy- 
droelastic analyses for which free surface strain energy and structural strain energy 
are equally significant. In most practical analyses involving liquid filled tanks the free 
surface strain energy is much smaller than the structural strain energy. Thus low 
frequency slosh dynamics is usually approximated with rigid structure, and flexible 
structure/fluid interaction is usually approximated with zero free surface pressure 
(the special case first developed in this section). 
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3 - HARMONIC REDUCTION OF GEOMETRICALLY AXISYMMETRIC STRUCTURES 

The NASTRAN hydroelastic analysis provides for a description of dynamics of 
axi symmetrically configured fluids in terms of circumferential harmonic pressure 
distributions. The distribution of pressure is 


+ 


N 

r 

k=l 



z. 

1 


P(r.,0., Zi) = Po(rj, z.) 
)cos z.)sin k0. 


(3.1-1) 


The fluid containing structure is described in terms of discrete physical dis- 
placements so that the structural representation is not limited to structurally axisym- 
metric containers. Coupling of harmonic pressure distributions with discrete struc- 
tural displacements in the NASTRAN formulation is not strictly consistent; moreover, 
in many cases it is inefficient. When a structure described in terms of discrete dis- 
placements is coupled with a fluid described in terms of circumferential harmonics, 
inconsistencies may arise if too few pressure harmonics are utilized; structural de- 
formation shapes associated with higher harmonics not included in the fluid represen- 
tation will reflect a lack of fluid inertia loading. Alternatively, when the discrete 
structural grid is too coarse to accurately describe the highest harmonic pressure 
distributions, large errors in the mode shapes associated with higher harmonics will 
be present. 


A consistent grid representation is realized by imposing a kinematic set of 
constraints on the structure restricting it to deform in circumferential harmonic pat- 
terns. For cylindrical coordinates (Fig. 3-la) the relationships between a discrete 
displacement and the harmonic displacement amplitudes are: 


Ui(ri,^i, z.) 



(3. l-2a) 
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Fig. 3-1 Cylindrical and Spherical Reference Frames 
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N 

V (r e ,z) = E 
^ ^ ^ ^ k=0 

U„ (r , z.)sin + U* (r., z.)cos k^ j 

_ Zk 1 1 X ^ ^ ‘J 

(3. l-2b) 

U3(ri,9 ,, Zj) = 

U„ (r , z )cos kd. + Uq (r , z )sin k0 
_ ^k 1 1 1 k ^ ^ \ 

(3. l-2c) 

N 

U(r 0 ,z)= E 
^ ^ ^ ^ k=0 


(3. l-2d) 

N 

U(r 0 z)= E 

^ ^ ^ ^ k=0 

[Ug (r., z.)cos k0. + U* (r., z.)sin k^^ 
1 k k 

(3.1-2e) 

N 

U(r 0 z)= E 

^ ^ ^ k=0 

fUg (r., z.)sin k0. + U* (r., z.)cos k0. 
Ik k 

(3. l-2f) 


For spherical coordinates (Fig. 3-lb) the relationships are: 


Ui (P ,e .,4>.) = E 

k ^ ^ ^ k-0 

U (P , 0 )cos k0. + U* (p,0.)sink<^ 

ik 1 1 1 1 1 1^ 

((3. l-3a) 

N 

u„ (P.,0 .,0.) = E 
\ ^ ^ ' k=0 

Uo (P.. 0.)cos k0 + U* (P.,0.)sink0 

Zk 1 1 1 Zj^ 1 1 

(3. l-3b) 

N 

u (P ,0 ,4>.) = E 
•^k ^ ^ ^ k-O 

fUg (P.,0 .)sin k0. + U* (P.,0 .)cos k^. 
1 k k 

(3. l-3c) 


j^U^k^P., 0 .)sin k<^. + U*k(/^.,0 .)cos k^. 

(3. l-3d) 
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U6^(p,,e,,^,) = ^r 


rUs^CPj.e j)sln k^j + U*^(P|,e j)cos k^J 

i)cos k^. + 9 psi" j 


(3. l-3e) 


(3. 1-3 f) 


Symmetric and antisymmetric displacement distributions are represented by un- 
starred and starred variables, respectively. It should be noted that all rigid body 
motions are represented by the above displacement distributions. 

Harmonic transformation of a discrete structural grid is accomplished in 
NASTRAN by the use of multipoint constraint (MFC) statements. The harmonic de- 
grees of freedom must be accounted for by supplementary grid points (in addition to 
the physical discrete grid). The general form of the harmonic transformation is 


where corresponds to the physical grid degrees of freedom to be transformed, 
and corresponds to the harmonic degrees of freedom (plus any discrete degrees of 
freedom not transformed). The transformation or constraint matrix Gra is composed 
of the appropriate sinusoidal functions evaluated at the discrete variable locations in 
accordance with the relationships outlined in Eq. 3. 1-2 and Eq. 3. 1-3. For a typical 
structure with JxK grid points such that there are J meridional rows and K circumfer- 
ential points in a row, the grid "g" set has typically 6xJxK degrees of freedom and the 
matrix semibandwidths are 6xK (assuming K< J). Application of the harmonic trans- 
formation as a reduction scheme where the number of harmonics, N, is much less than 
K results in a Uj^ set of 6xJxN generalized coordinates with matrix semibandwidths of : 
6xN. If N«K harmonic reduction represents a radical reduction in the number of de- 
grees of freedom as well as matrix bandwidth. Further reduction of the system de- 
scription is possible by a small Guyan reduction by choosing the generalized rotation 
degrees of freedom and tangential displacements as members of the omitted set of dis- 
placements. In such a case the analysis set consists of JxN degrees of freedom. This 
represents a radical reduction in degrees of freedom by a factor of (N/6K) without a 
costly large matrix decomposition typical of Guyan reduction. 
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The economy and accuracy of harmonic reduction was first tested on the spheri- 
cal cap with uniform thickness illustrated in Fig. 3-2. A nodal grid consisting of 20 
circumferential divisions in a semicircle and 10 meridional divisions was chosen 
resulting in a 1266 DOF "g-set" structural model. Three circumferential harmonics 
(0, 1, 2) were chosen for harmonic reduction. It should be noted that the apex node is 
left in terms of rectangular coordinates since the p>olar degrees of freedom have no 
meaning at this node. After application of the fixed-base boundary condition and sj^m- 
metric kinematic constraints at the pole the reduced representation consists of a 138 
DOF "f-set". A small Guyan reduction omitting the non-zero displacements at the 
pole and rotational and circumferential generalized displacements yields a 72 DOF 
"a-set". At this point, all natural frequencies and the first 15 modes were calculated 
by the Givens method. 

The results of the above strategy were then compared to STARS-II (Ref. 8) 
results which were assumed exact. In addition, NASTRAN results utilizing various 
Guyan reduction strategies illustrated in Figs. 3-3 and 3-4 were compared to the 
STARS-II and harmonic reduction results. 

A comparison of natural frequency results (Table 3-1) indicates that the overall 
accuracy of the 72 DOF harmonic reduction representation is better than the 190 DOF 
Guyan reduction representation. The quantities are representative of fractional 
frequency errors and A 2 , the frequency squared errors, are representative of mode 
shape errors. Comparisons of the first two axisymmetric mode shapes illustrated 
in Fig. 3-5 Indicate that although frequencies are rather accurate, mode shapes may 
contain large local errors. Mode shapes resulting from harmonic reduction are 
clearly more accurate than those resulting from Guyan reduction. 
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Fig, 3-2 60^ Spherical Cap Model, (h/R) = 0.05 
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Fig. 3-3 Analysis Set Nodes for Guyan Reduction G-110 
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Fig. 3-4 Analysis Set Nodes for Guyan Reduction G — 95 and G — 190 


Table 3-1 60® Spherical Cap *- Comparison of Modal Frequencies 

Aj = , A 2 = (co/a?g)^ — 1 

STARS Harm. Red. H 72 Guyan Red G110 GUYAN Red. G -95 GUYAN Red. G 1^ 
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Fig. 3-5 60° Spherical Cap, Fixed Base Mode Shape Comparisons 
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Another aspect of the harmonic reduction technique is its relative efficiency. A 
comparision of harmonic and Guyan reduction Central Processor Unit (CPU) times pre 
sented in Table 3-2 indicates a significantly lower reduction CPU time for harmonic 
reduction. This is attributed to elimination of a large scale matrix decomposition, 
characteristic of the Guyan reduction, and to the fact that fewer degrees of freedom 
are required for comparable accuracy (e. g. , harmonic 72 DOF-5% accuracy vs. Guyan 
190 DOF-12% accuracy). The CPU time associated with eigenvalue analysis of the har- 
monic analysis set is naturally much less than that associated with the Guyan reduction 
analysis set. 

For the current structural model, circumferential harmonics were uncoupled 
due to the axisymmetry of shell thickness. In cases where thickness varies with cir- 
cumferential location ( 0 ) harmonics will be coupled. The degree of coupling is a func- 
tion of the relative abruptness or smoothness of thickness distribution and the thick- 
ness asymmetry. 

A general purpose FORTRAN IV computer program for generation of multipoint 
constraint bulk data cards, HARM, has since been written for typically large hydro- 
elastic problems, many constraint statements must be prepared. A listing and de- 
scription of this program are presented in Appendix A of Volume II. 

Table 3-2 60° Spherical Cap Modal Accuracy/Cost Comparisons (1266 DOF Grid Set) 


Reduction 

Scheme 

Analysis Set Size 

Modal Accuracy* 
(15 Modes) 

Total Central Processor 
Unit Time (Sec) 

Harmonic 

72 

0.053 

238 

Guyan 

110 

0.276 

346 

Guyan 

95 

0.149 

328 

Guyan 

190 

0.117 

531 

N 

' r, .o 12 



’Error = -j!j 1 




i 

= iL^ j 
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4 - NUMERICAL RESULTS 

The present hydroelastic analysis has been implemented in NASTRAN and veri- 
fied and demonstrated on a number of problems. The problems fall into two catego- 
ries, namely, analytical verification problems for which exact solutions are known and 
demonstration problems for which experimental data are available. The l/8-scale 
shuttle external tank is included in the second category. 

4. 1 ANALYTICAL VERIFICATION PROBLEMS 

Ex. 1. Fluid Filled Hemispherical Container 

This first problem consists of the fluid filled hemispherical container 
illustrated in Fig. 4-1. The container is massless and follows the artificial structural 
law 


p = au (4.1-1) 

r 

The exact free vibration solution for this problem is expressed in terms of spherical 
harmonics (Ref. 9). The modal displacements on the structural surface and free sur- 
face are 


U = P^(cos 0 )cos n cp 


(4. l-2a) 


, „ dP” (cos 0 ) 

%=lRj 


(4. l-2b) 


respectively. The modal pressure function is 


P = 




mn 


m 


m 

P^(cos 0 )cos n(p 


(4. 1-3) 


and the natural frequencies are 
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Fig. 4-2 Fluid in a Hemispherical Container — Finite Element Model 
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The function P^(cos 0 ) is the associated Legendre function for which m is the merid- 
ional wave index and n is the circumferential wave index. The allowable indices are 
such that the sum m+n must be odd with m > n. 

The finite element model of the hemispherical fluid filled container is illustrated 
in Fig. 4-2. As a result of the structural law Eq. 4. 1-1 the structural stiffness ma- 
trix is diagonal with entries 

K.. = OJA. (4.1-5) 

11 1 

where A. is the area associated with the "ith" (radial) degree of freedom. All diago- 
nal entries not associated with radial degrees of freedom are null. The fluid model is 
expressed in terms of the circumferential pressure harmonics n = 0, 2, 4 and the 
structural surface and free surface grids are reduced by harmonic reduction accord- 
ingly. The fluid mass matrix developed in the modified Rigid Format 7 version is ex- 
pressed in terms of a 21-degree of freedom analysis set of structural radial displace- 
ments (7 meridional locations, circumferential harmonics n = 0, 2, 4). 

Natural frequencies and mode shapes for the finite element model were calcu- 
lated in Rigid Format 3 by the Givens method. A comparison of exact and NASTRAN 
calculated nondimensional natural frequencies is presented in Table 4-1 and compari- 
sons of selected modal displacement distributions are presented in Fig. 4-3. In gen- 
eral the finite element results are in excellent agreement with the exact solution; as 
expected in any finite element analysis the level of accuracy decreases with modal 
complexity. 

Ex. 2. Fluid Filled Circular Cylindrical Shell 

The second verification problem consists of the fluid filled circular cy- 
lindrical shell illustrated in Fig. 4-4. The shell structure is taken as one with bending 
as well as membrane stiffness. The geometric properties of the shell consist of a cy- 
linder aspect ratio Z/R = 2 and a thickness ratio h/R = 0. 01. In addition, the fluid to 


4-3 



Volume I 


structure density ratio is Pf/Ps = 1/3 and the structural material Poisson ratio is 
= .3 An exact hydroelastic modal solution is known for an infinitely long cylinder 
which holds for the present problem when the structure is subjected to the boundary 
conditions 


Ur = = 0 (shear diaphragm) @ z = Z, r = R (4. l-6a) 

P = 0 (free surface) @ z = Z, r<R (4. i-6b) 

dz = N = 0 (symmetry) @ z=0, r=R (4. l-6c) 

z " 

U^ = 0 (fixed bottom) @ z=0, r<R (4. i-6d) 

The exact free vibration solution is expressed in terms of cylindrical harmonics 
(Ref. 10). The normalized modal displacements on the structural surface (r = R, 

0< z<Z) and the free surface (r<R, z = Z) are 


U 

z 


m 7TZ 

cos 2'z cos ne (4. l-7a) 

_ mTT j mTTr 
^ 2Z n 2Z 



cos n 0 

r=R 


(4. l-7b) 


respectively, where (^f^) ^ modified Bessel function. The modal pressure 

function on the structural surface is 



(4. 1-8) 
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Table 4-1 Fluid in a Hemispherical Container Natural Frequency Comparisons 
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Fig. 4-3 Fluid in Hemispherical Container - Mode Shape Comparisons 
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Noting that the structural deformation shape Eq. 4. l-7a is the same as for an empty 
cylinder the hydroelastic natural frequencies are related to the empty structural fre- 
quencies fjv as 
mn 

“e 

mn 

(4. l-9a) 




where mf and mg represent the effective fluid and structural masses 

^f n\ 2Z / 

~ f ^ /t /mTrr \ \ 
\2z/^[n\^2Z 


nig = Psh ( 4 . i-9c) 

The empty shell frequencies are known exactly on the basis of Flugge's shell 

theory (Ref. 11). 

The finite element models of the shell and fluid are illustrated in Fig, 4-5, The 
structural grid for the quarter shell consists of 924 degrees of freedom and the fluid 
grid consists of 165 degrees of freedom (55 nodes of rotation, circumferential har- 
monics n = 0, 2, 4). The harmonic transformation retaining harmonics n = 0, 2, 4, 
application of single point constraints to enforce boundary conditions, and a small 
Guyan reduction retaining only radial displacements ultimately resulted in a 30-degree 
of freedom analysis set. Listings of the NASTRAN Rigid Format 7 and Rigid Format 
3 data for this problem are presented in Appendix B of Volume II. 

All natural frequencies and 25 mode shapes with and without th6 fluid included 
were calculated in Rigid Format 3 by the Givens method. The mode shapes in both 
cases were nearly identical to one another as concluded in the exact analysis. Selected 
mode shapes for the liquid filled case are illustrated in Fig. 4-6. Frequency spectra 
for the empty and fluid filled shells are presented in Fig. 4-7 illustrating excellent 
comparison between finite element NASTRAN and exact results in both cases. 
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= 30 DOF 


Fig. 4-5 Circular Cylindrical Shell with Fluid — Finite Element Model 
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Fig. 4-7 Circular Cylindrical Shell Frequency Spectra 
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A characteristic of the modified NASTRAN hydroelastic analysis which is as 
significant as numerical accuracy is computational economy. On the Grumman IBM 
370/165 computer the Rigid Format 3 solution time for the empty cylinder was 2 min. 
2 sec; for the fluid filled cylinder the Rigid Format 7 and 3 solution times were 1 min 
15 sec and 2 min 6 sec, respectively. These computation times represent consider- 
able cost savings for NASTRAN hydroelastic analysis and empty-structure modal 
analysis. 

Ex. 3. Pressurized Fluid Filled Circular Cylindrical Shell 

The above described fluid filled cylindrical shell is now considered in 
a uniformly pressurized state. A closed form free vibration solution in this case 
follows Eq. 4. 1-6 through 4. 1-9 withWg^^ representing the empty pressurized shell 
frequencies. An approximate expression for the adjusted empty shell frequency with 
pressurization, P^, (Ref. 10) is 


= 1 + • 


u 


E 


Eh 


mn 



(4. 1-10) 


The mode shapes in the pressurized case are the same as in the unpressurized 

case. 


A modified version of Rigid Format 13 (normal modes with differential stiffness) 
was used to calculate modal data for the empty and fluid filled cylinder representation, 
with a pressurization level (1— v )PQR/Eh = 0. 001. The models used for this analysis 
and all reductions are the same as for the unpressurized case. The frequency spectra 
for the pressurized empty and fluid filled shells are presented in Fig. 4-8 for the n=4 
modes only; pressurization caused negligible changes in the n = 0, 2 modes. The na- 
tural frequencies resulting from the NASTRAN analysis are in excellent agreement 
with the results based on Eq. 4. 1-10. In addition, as in the unpressurized analysis, 
computation times were quite satisfactory. Rigid Format 13 modal analysis CPU 
times were 3 min 2 sec and 3 min 20 sec for the empty and fluid filled cylinders, re- 
spectively; approximately 1 min CPU time was required to generate the differential 
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stiffness terms in Rigid Format 13. Rigid Format 13 data for this problem are pre- 
sented in Appendix B of Volume II. 

4. 2 COMPARISONS WITH EXPERIMENTAL DATA 

Ex. 1 Liquid Filled Cylinders Under Static Pressurization 

A detailed experimental study of the dynamics of structurally axisym- 
metric and asymmetric circular cylinders under various liquid (water) fill and static 
pressurization conditions has been conducted at NASA Langley Research Center by 
Mr. Robert Herr. Data resulting from these tests (unpublished) are very complete 
and serve as excellent information for analysis/test correlation studies. The test 
articles are aluminum cylinders with mean radius 25.4 cm and height 50. 8 cm. The 
cylinder walls are welded at the top and bottom to heavy aluminum plates as illus- 
trated in Fig. 4-9. The axisymmetric test article has a cylinder wall thickness of 
. 08128 cm and the asymmetric test article has a wall thickness variation around the 
circumference 0. 0508 — 0. 1016 according to the equation 


h 

h 

max 


= 0. 75-0.25 cos 


e 


(4.2-1) 


The NASTRAN finite element fluid and structural grid representations for the 
1/2 fill condition are illustrated in Fig. 4-10. The structural model for a l/2 cylinder 
(0°< 6 <180°) is described by a sufficiently fine grid to simulate the higher circumfer- 
ential harmonic shapes (e.g. n = 0-15) which are known a priori from the experimental 
results to dominate in the lowest frequency modes. The fluid representation consists 
of a fine radial grid near the structural wall to insure simulation of the sharp pressure 
gradients in the higher circumferential harmonics; the generalized fluid pressures are 
expressed in terms of the symmetric harmonics (cos n^) for n — 0 to 15. The grid set 
for the half filled cylindei- consists of 480 pressure degrees of freedom and 2,046 struc- 
tural degrees of freedom. The structural description in this case is not reduced by 
harmonic reduction since all harmonics up to n = 15 are of interest and insignificant 
computational economization would be gained by the harmonic transformation. 


4-13 





> 



QC 


> 

O 


X 

UJ 


o 

I 


UJ 

I- 


I 

Q 


h- 

Ul 


Q 

UJ 

N 

X 

N 

q: 

D 

'{k 

Z 

o 

1 “ 

13 

<o 

CO 

UJ 

QC 

-J 

O 

UJ 

Q. 

CO 

QC 

z 

z 

0 . 

D 

< 

H 

I— 

QC 

O 

O 

I— 

< 

< 

CO 

X 

X 

< 

UJ 

UJ 

z 

CO 

CO 

CO 

UJ 

UJ 

UJ 

H 

!- 

H 

O 

O 

O 

Z 

z 

z 

UJ 

UJ 

UJ 

a 

Q 

Q 

1 

1 

• 


(Q 

h. 

ts 

K 

00 

t? 

c 

V 

3 

O’ 

0 > 

w 

U. 

II 

c 

£ 

00 


■o 

> 

o 


■a 

V 

N 


0 ? 

w 

LL 


4-14 






■ ^0.0813 CM LIP (TOP & BOTTOM) 



50.8 CM 




50.8 CM DIAMETER 


Fig. 4-9 Circular Cylindrical Shell Test Article 
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The first series of problems studied pertained to the cylinder of uniform thick- 
ness. The empty cylinder was first considered with an assumed axial plane of symme- 
try at z=25.4 cm such that only m=l, 3, 5 modes would be calculated. The grid set of 
1116 degrees of freedom consisting of nodes below z=25.4 cm was reduced by Guyan 
reduction to an analysis set of 276 radial degrees of freedom with the base assumed 
completely fixed (clamped). The analysis on Rigid Format 13 was performed for the 
test conditions of zero pressurization and static pressurization at 5.516xlO‘^N/m2. The 
frequency spectrum of m = l, n>4 modes is illustrated in Fig. 4-11 along with the test 
results. The calculated frequency spectrum was higher in both cases than the experi- 
mental frequency spectrum. A series of structural model modifications to reconcile 
the differences in results were considered and it was finally concluded that axial flexi- 
bility in the cylinder/plate weld provided the proper correction. Incorporation of the 
boundary conditions 

U = = N = 0 @ z = 0, 50. 8 cm (4. 2-2) 

r ^ z z 

4/2 

resulted in extremely accurate frequency spectra for the empty 0 and 5. 516x10 N/M 
pressurization conditions, respectively, as illustrated in Fig. 4-11. 

The half filled, unpressurized condition was then considered. The liquid free 
surface was described in terms of single point constraints applied to the surface pres- 
sures; free surface displacements were not desired as output information. Retaining 
only the radial displacements below the free surface in the analysis set a 248 degree 
of freedom fluid mass matrix and pressure recovery matrix were calculated with the 
modified version of Rigid Format 7. The cylinder structure in this case does not have 
a dynamic plane of symmetry at z = 25. 4 cm; the lower portion (z<25. 4 cm) is loaded 
by the fluid inertia and small structural inertia while the upper portion (z >25. 4 cm) is 
loaded only by the structural inertia. This provides motivation for Guyan reduction 
with all degrees of freedom at and above z = 25. 4 cm (not including the supported de- 
grees of freedom) "omitted". A Guyan reduction on the structure was then performed 
resulting in an analysis set consisting of 248 radial degrees of freedom. Hydroelastic 
modes based on the clamped and modified clamped Eq. 4. 2-2 end conditions were then 
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calculated. The m = l, n>4 modal frequency spectra illustrated in Fig. 4-11 show 
that the representation with the modified end conditions is quite accurate as concluded 
in the empty cases. At this point it was felt that a hydroelastic analysis with pres- 
surization was not required since enough confidence in the model was gained from the 
three cases considered above. A study of the unsymmetric cylinder dynamics was then 
initiated. 

The unsymmetric cylinder structural model consists of the same grid set as in 
the case of the axisymmetric cylinder with circumferential thickness variation. The 
hydroelastic study of this cylinder was limited to the half filled condition with 0 and 
5. 516x10'^ N/m^ pressurization and the realistic edge condition Eq. 4.2-2 was ap- 
plied. The m = l mode shapes calculated in Rigid Format 13 for the unpressurized and 
pressurized conditions illustrated in Figs. 4-12 and 4-13, respectively, are in excel- 
lent agreement with the test results as are the modal frequencies (presented in the il- 
lustrations). 

Computation times for the cylinder study were moderate since harmonic reduc- 
tion was not appropriate and thus not utilized. In all cases considered, all eigenvalues 
and 25 eigenvectors were calculated by the Givens method. Computation times for the 
empty axisymmetric cylinder (1116 DOF g-set, 276 DOF a-set) were 509 CPU sec and 
524 CPU sec for the unpressurized and pressurized cases, respectively. Preparation 
of fluid matrix data in Rigid Format 7 required 97 CPU sec and computation of hydro- 
elastic modes (2,046 DOF g-set, 248 DOF a-set) required 1,193 CPU sec and 1254 CPU 
sec for the unpressurized and pressurized cases, respectively. The increased CPU 
time required is predominantly due to the increased structural grid set size of the fluid 
filled cases; the increase in Guyan reduction time for systems of equivalent matrix 
bandwidth is proportional to the increase in g-set degrees of freedom. Computation 
times for the unsymmetric cylinder were similar to those required for the axisymme- 
tric cylinder. 

Ex. 2 l/8-Scale Space Shuttle External Tank 

An investigation of 1/8 scale space shuttle external tank dynamics in a 
free-free supported condition has been undertaken. The 1/8 scale external tank 
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consists of two separate propellant tanks connected by a cylindrical section. The finite 
element hydroelastic model described in detail in Volume II and Ref. 12 consists of a 
grid set of 348 pressure degrees of freedom, 2,058 structural degrees of freedom and 
768 harmonic structural degrees of freedom. Harmonics n=0, 1,2,3 were chosen to 
describe asymmetric dynamics with the pitch plane taken as an axis of symmetry. 

The analysis set of displacements resulting from a combination of harmonic and 
Guyan reductions consists of 128 harmonic degrees of freedom associated with out- 
ward normal motion of the tank wall. 

Three liquid fill conditions have been studied consisting of liftoff, post max Q 
and empty. In terms of liquid height above the respective bulkheads, the conditions 
are identified as: 

• liftoff h = 190. 5 cm, hj^jj = 358.14 cm 

2 

• post max Q: ^LH =330.2 cm 

2 

The liquid free surfaces are taken normal to the tank axis thus ignoring some 
free surface tilt to be experienced in flight since the NASTRAN hydroelastic analysis 
is currently limited to axisymmetric fluid configurations. This limitation, however, 
can be overcome by utilization of heat conduction polyhedral finite elements as fluid 
elements according to the analogy presented in Appendix B. 

For each of the fill conditions, 128 natural frequencies and 25 mode shapes and 
modal pressure distributions were calculated with very good computational efficiency. 
About 20 CPU minutes per liquid level on the Grumman IBM 370/165 computer was 
required to perform the entire analysis including matrix assembly, reduction and 
modal analysis. In previous attempts to study the dynamics of the same finite element 
representation with the old NASTRAN hydroelastic analysis, computation times were 
In excess of 70 CPU minutes with only one natural frequency and mode shape computed. 
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Modal data for each of the three fill conditions studied are summarized in 
Tables 4-2 through 4-4 with dome pressure gain presented as a measure of relative 
significance in Pogo susceptibility. The pressure gain is defined as the dome pres- 
sure (LOg and LHg) in the vicinity of the feedline interface near the respective tank 
bottoms resulting from a unit modal acceleration. Plots of the mode shapes corre- 
sponding to the intermediate post max Q fill condition are presented in Figs. 4-14 
through 4-35. 

Excellent agreement between analysis and experimental frequencies occurred 
in the first axial mode but poor agreement occurred in the bending modes. The 
source of the discrepency is believed to be in the finite element representation of the 
structure and efforts to resolve the discrepency are discussed in Volume II. 


Table 4-2 1 /8-Scale External Tank Hydroelastic Mode Summary (at Liftoff) 


Mode No. 

Freq. (Hz) 

Modal Mass 

Description of Mode 

LO 2 Dome 
Pressure Gain 
X 103 

LH 2 Dome 
Pressure Gain 
X 103 

4' 

29.7 

4,751 

ET 1st Axial n=0 

0.130 

0.034 

5 

34.5 

0.857 

LO 2 n=2 (No Dome) 

0.003 

0.001 

6* 

35.7 

0.760 

ET 1st Bending n=1 

0.040 

0.017 

7 

36.6 

0.428 

LO 2 n=3 (No Dome) 

0.008 

0.005 

8' 

54.9 

2.667 

ET 2nd Axial n=0 

0.067 

0.091 

9' 

57.8 

0.131 

LH 2 Cylinder n=2,3 

0.016 

0,019 

10’ 

61.4 

0.067 

LH 2 Cylinder n=3,2 

0.017 

0.034 

11 

62.1 

0.395 

LO 2 n=3 (No Dome) 

0.011 

0.005 

12* 

63.8 

0.520 

ET 2nd Bending n=l 

0.070 

0.006 

13 

68.4 

0.581 

LO 2 n=2 (No Dome) 

0.002 

0.003 

14* 

96.0 

0.618 

LO 2 n= 1 

0,102 

0,001 

15* 

96.1 

0.433 

LO 2 n^O 

0.302 

0,006 

16* 

109.4 

0.455 

LO 2 , LH2n-0 

0.155 

0.028 

17 

114.5 

0,741 

LO 2 , LH2n=2,3 

0.009 

0.001 

18 

117.8 

0,277 

LO 2 n=3 (No Dome) 

0.006 

0 

19* 

119.7 

0,142 

LH 2 Cylinder, LOX n-2,0 

0.085 

0.002 

20* 

119.8 

0.254 

LO 2 n-0 

0.276 

0,017 

21* 

124.2 

0.221 

LO 2 n=l 

0.253 

0,003 

22 

128.6 

0.062 

LH 2 Cylinder n=3 

0.002 

0 001 

23* 

135.0 

0.475 

LO 2 n=0 

0.150 

0.025 

24* 

135.9 

0.431 

ET. LOX Dome n=1,0 

0.194 

0,005 

25 

138.2 

0.589 

LO 2 n=2 

0.006 

0,002 


'Denotes POGO Sensitive Mode 

NOTE: Modes 1, 2, 3 are Rigid Body Pitch Plane Modes 
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Table 4-3 1 /8-Scale External Tank Hydroelastic Mode Summary (Post Max Q) 
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Table 4-4 Empty 1/8-Scale External Tank Mode Summary 


Mode No. 

Freq. (Hz) 

Modal Mass 

Description of Mode 

4 

105.6 

0.0425 

ET 1st Bending n=1 

5 

153.0 

0.0094 

LH2 Cylinder n=3 

6 

161.7 

0.0172 

LH2 Cylinder n = 2, 3 

7 

226.0 

0.0497 

ET 2nd Bending n=1 

8 

257.8 

0.0770 

ET 1st Axial n=0 

9 

274.7 

0.0271 

LH2 Cylinder n=2,3 

10 

328.3 

0.0122 

LH 2 Cylinder, LOX n=3,2 

11 

332.0 

0.0149 

LO 2 , LH2 n=3,2 

12 

332.8 

0.0234 

LO 2 , LH 2 n=3,2 

13 

343.7 

0.0118 

ETn=3 

14 

357.8 

0.0696 

ET Bending n=1,3 

15 

431.0 

0.0210 

LH 2 Cylinder n=2 

16 

459.1 

0.0615 

LH 2 Cylinder n=3, ET n=l 

17 

472.9 

0.0114 

LH 2 Cylinder n=3 

18 

482.2 

0.0185 

ETn=2 

19 

498.6 

0.0076 

LO 2 n=3 

20 

513.2 

0.0697 

ETn=l,2, 3 

21 

533.1 

0.0243 

ET n=2, 1 

22 

567.2 

0.0487 

ETn=1,2,3 

23 

604.6 

0.0391 

ETn=2, 1,3 

24 

625.4 

0.9116 

LO 2 n-3 

25 

628.0 

0.0144 

ET n=3,2 


4-25 














Volume I 



4-26 


Fig. 4-14 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 4 


Volume 1 



4-27 


Fig, 4-15 1/8-Scale Space Shuttle External Tank — Post Max Q, Mode 5 
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Fig, 4-16 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 6 
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Fig, 4-17 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 7 


Volume I 



4-30 


Fig. 4‘18 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 8 
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Fig. 4‘19 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 9 
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Fig. 4-20 l/8-Scale Space Shuttle External Tank - Post Max Q, Mode 10 
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Fig. 4-21 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 11 
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Fig. 4-22 1/8-Scale Space Shuttle External Tank — Post Max Q, Mode 12 
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Fig. 4<23 1/8-Scale Space Shuttle External Tank — Post Max Q, Mode 13 
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Fig. 4-24 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 14 
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Fig. 4-25 1/8-Scate Space Shuttle External Tank - Post Max Q, Mode 15 
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Fig. 4-26 1/8^cale Space Shuttle External Tank — Post Max Q, Mode 16 
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Fig. 4-27 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 17 
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Fig. 4-28 1/8-Scale Space Shuttle External Tank — Post Max Q, Mode 18 
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Fig. 4<29 1/8-Scale Space Shuttle External Tank — Post Max Q, Mode 19 
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Fig, 4-30 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 20 
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Fig. 4-31 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 21 
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Fig. 4-32 1 /8-Scale Space Shuttle External Tank — Post Max Q,Mode 22 
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Fig. 4^33 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 23 
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Fig. 4*34 1/8*Scale Space Shuttle External Tank — Post Max Q, Mode 24 
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Fig. 4-35 1 /8-Scale Space Shuttle External Tank — Post Max Q, Mode 25 



Volume I 


5 - CONCLUSIONS 

Symmetric finite element matrix formulations for compressible and incom- 
pressible hydroelasticity have been developed on the basis of Toupin's complementary 
formulation of classical mechanics. The incompressible formulation applicable in 
propellant tank hydroelastic analysis has been implemented in NASTRAN to replace 
the unsymmetric matrix formulation. The new technique which utilizes existing fluid 
and structural finite elements has been verified and demonstrated to be accurate and 
efficient. 

The fluid representation according to the new technique consists of a symmetric 
fluid mass matrix described in terms of surface deformation only and an additional 
surface stiffness matrix when gravitational potential and ullage pressure stiffness are 
Included in the fluid idealization. The fluid mass and stiffness matrices are then 
added directly to the structural mass and stiffness matrices, respectively, forming a 
symmetric set of hydroelastic equations in terms of structural displacements. As a 
result of the extensive NASTRAN structural modeling capability, differential stiffness 
effects due to static pressurization and fluid weight may be accounted for in the struc- 
tural idealization. Modal hydroelastic analysis is performed with the same efficiency 
as in the case of a non-fluid filled structure as a result of very few additional degrees 
of freedom being required for the fluid. 

The efficiency of the new hydroelastic analysis technique has been enhanced for 
both fluid and structure by introduction of harmonic reduction, applicable to geomet- 
rically axisymmetric structures, as an alternative to Guyan reduction. When the 
number of harmonics utilized is much less than the number of discrete nodes about a 
circumference, overall matrix size and bandwidth are significantly reduced. Har- 
monic reduction which was developed early in the present study was first demonstrated 
on a 60° spherical cap exhibiting superior accuracy and efficiency relative to Guyan 
reduction for an axisymmetric structure. 

The formulation has been verified by comparison with exact analytical results 
for a fluid filled hemispherical container, a fluid filled circular cylindrical shell and 
a pressurized fluid circular cylindrical shell. In all three cases, excellent correla- 
tion was exhibited as well as very good computational efficiency. 
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Analysis/test discrepancies on the 1/8-scale external tank model have not yet 
been resolved. The efficiency of the current 1/8-scale external tank analysis, 
however, is very encouraging upon comparison of computation times between the 
present analysis and the standard unsymmetric NASTRAN hydroelastic formulation. 
Approximately 20 CPU minutes of computer time is required to extract 128 natural 
frequencies and 25 mode shapes with the present anal 3 d;ical technique while more 
than 70 CPU minutes was required with the standard NASTRAN hydroelastic analysis 
to extract a single natural frequency and mode shape. 

The analysis/test correlation study on symmetric and imsymmetric circular 
cylindrical shells under various fluid fill conditions and static pressurization is 
considered very good. 

It is strongly recommended that the present NASTRAN hydroelastic analysis be 
extended to include the capability for modeling non-axlsymmetric fluid geometries by 
use of polyhedral heat condition finite elements. Long term goals anticipating future 
analytical requirements should also be pursued including Implementation of the 
symmetric compressible hydroelastic formulation and development of general 
polyhedral compressible fluid elements. 
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APPENDIX A - TOUPIN'S VARIATIONAL PRINCIPLE 

A complementary variational formulation of classical mechanics utilizing 
impulse quantities as generalized coordinates was introduced by Toupin in 1952 
(Ref. 6). This formulation, which may be interpreted as a dynamic generalization of 
the complementary energy principle of statics, has not received much attention in the 
development of analytical methodology; it is advantageous in most analyses to utilize 
displacement quantities as generalized coordinates. As in the displacement formula- 
tion of mechanics, complementary counterparts of D'Alembert's Principle, Hamilton's 
Principle and the Euler- Lagrange equations are derivable as consequences of Newton's 
laws. The development presented here closely parallels the derivation of the varia- 
tional displacement formulation in mechanics. 

Newton's second law states that for a particle 

4-(m“r)=fi^ (A-1) 

with^ and F representing the velocity and force vectors, respectively, and with m^ 
and taken as constants. Newton arbitrarily defined m^ as the particle mass and 
as unity thus setting the course for development of displacement oriented formula- 
tions in mechanics. By arbitrarily defining m^ as unity and as the particle inert- 
ance the course is now set for development of alternative formulations; the inertance 
is simply the Inverse of mass. Integration of (A-1) with respect to time results in 

t 

T(t) = I ^ of ^dt + r(o) (A- 2a) 

o 

or by taking the initial velocity as a consequence of all previously applied force 

o 

r{o) = J F dt (A-2b) 

-00 

the concise integrated statement of Newd;on's law is 

T(t)=f^F(t) (A- 2c) 
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with ^ 

F(t) = j F(t) dt (A- 2d) 

- eo 
A 

The quantity F(t) is the total impulse which has brought the particle from rest to the 
current velocity. The above statements also hold for aggregates of particles. 

It is now postulated that a system of n particles moves such that under the 

A 

arbitrary virtual impulses, 5F., 


£ <^i- 'Vo (A- 3 ) 

i=l 

where r. denotes the position vector with respect to a Newtonian reference frame. 
The above is the complementary counterpart of D'Alembert's principle. Integration 
of the above expression over the time interval (t^, tj^) results in 


1 n A ;i n A A 

I Eir.-6F.)6t- I E F. • <}F.dt=0 

t i=l t i-1 

o o 


(A- 4) 


At this point it is convenient to define the complementary kinetic energy and virtual 
work fimctions which are 


n 


A A 


T =1/2 E rt . F. • F. 
c i^oi 1 1 

i=l 


(A-5a) 


and 


n 


— A 


6W = E r' 6F. 
c i 1 

i=l 


(A-5b) 


respectively. The time integral of the virtual work function integrated by parts is 


r / n ^ A 

I dW^dt= r E (r. •<5F.)dt 


t i=l 
o 
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n 

I 

1=1 


A 


r. 

I 


6^; 


1 n _U A 

r 2 : (^ 1 * 


By imposing the constraints 

dF. (t^) = 6F. (t^) = 0 


(A-6a) 


(A-6b) 


the complementary virtual work function is expressible in the alternate form 


5 W = - Y (^ • d F.) (A-6c) 

C A A 

i=l 

Incorporation of Eq. A-5a and Eq. A-6c into Eq. A-4 yields the statement of the 
complementary Hamilton's principle 

/i 

/ (dT + dW ) dt = 0 (A-7) 

J 0 c 

t 

o 

If the complementary work function, W , is expressible as 

c 

w = - U + w/ (A-8) 

c c c 


where U defined as the complementary potential energy function has the general 
c 

fimctional dependence 


U = U 
c c 


A 

(F, 


A A 

Fn» F^ »■ 


A 

, F 


n 


;t) 


(A- 9) 


and W"^is the remaining part of the work function best thought of as dependent on 
c 

externally applied displacements (or velocities) such that 


(5W^=- X (T? • dF.) (A-10) 

c ^ 

i=l ® 

The complementary Lagrangian function is now designated as 

L = T - U (A-12) 

c c c 
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and the complementary form of Hamilton's principle is 

A A A A 

sj Fn, t)dt+ J <jw^dt = 0 (A-12) 


In order to derive a general set of complementary Euler- Lagrange equations 
(Toupin's equations) the notation of generalized impulse coordinates is introduced. 
The generalized impulses are defined as those which 

• Determine the dynamic configuration of the system 

• May be varied arbitrarily and independently without violating the constraints 
of the system. 

Consider now a system of n particles for which the individual impulses are fvinctions 
of m (m < n) generalized impulses, i.e.. 


A A A A 

F. = F. (Q^ , . . . , Q^, t), i = 1, , 


, n 


(A-13) 


A 

The functions F. are assumed continuous with continuous partial derivatives and the 

1 A 

time derivative of the total impulse F^ is 


A A A 

A dF A 9F a 9F 

Fi= a^Qi+--*+ — + 

5Q, dQ ™ 

1 m 


9t 


Thus the complementary Lagrangian has the fimctional dependence 


A A A A 

L — L (Q- , . . . , Q , Q- ,...,Q »t) 

c c '1 ’ m’ 1 ’ m’ ' 


(A-14) 


(A-15) 


A 


The variation of impulse, F^ is 


A A 

,A A ,A 


m 


(A- 16) 


m 


and the external complementary virtual work Eq. A- 10 is expressed as 


dW''= - rq- Q, +.. .+ q 1 

c Ig 1 m^ m J 


(A- 17a) 
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where 



Substitution of Eq. A-15 and Eq. A-17 into Eq. A-12 results in 



(A-17b) 


(A-18) 


Integration by parts noting the constraints at the limits 

(y = 6Q. (t^) = 0 


yields 



(A-19) 


(A- 20a) 


where the Euler- Lar^range equations are the integrands 


dt 


8L 


a 


dL 

7^ + q 


a 


je 


= 0 


(A- 20b) 


Application of the complementary principle to large deformation and other non- 
linear problems requires a re-education of the analyst; one's physical intuition must 
be reoriented towards kinetic rather than kinematic considerations. Restricting 
further discussion to small deformation problems (in particular conservative small 
vibration problems), the kinetic and potential energy ejqjressions reduce to the 
quadratic forms 


T = T 
c c 




) 


2 


m 

E 

1=1 


m 

E 

j=l 


A A 

L. . Q. Q, 

ij ^1 


(A- 21a) 
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mm A 

= 2 £ X 
i=l j=l 


(A-21b) 


respectively. The matrix set of Euler- Langrange equations resulting from Eq. A-22 


(L) |q| + (C) = |q^J 


Taking the time derivative (noting that Q = Q) the result is 


(L) <Q^ + (C) 


f'^} = {'•e} 


(A- 22a) 


(A- 22b) 


where |q| is the generalized force vector, (L) is the generalized inertance matrix 
and (C) is the generalized flexibility matrix. The significance of the matrices becomes 
apparent when the equations are rearranged as 


where 


(L) {q} = + {q^} = } 


(A- 23) 


(A- 24) 


are interpreted as "complementary D'Alembert accelerations". Thus (L) is the 
inverse of the generalized mass matrix (M). In the case of statics the kinetic energy 
is null and the special equation (integrated twice) 


H = K)> 


(A-25) 


represents a statement of static equilibrium; and the flexibility matrix (C) is identified 
as the inverse of the generalized stiffness matrix. 
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The qualities of a complementary formulation are best illustrated by a set of 
simple examples. 

Example 1: 

Consider first the simple oscillator with base excitation, U^, illustrated in 
Fig. A-1. The kinetic and strain energies are 

respectively and the complementary virtual work is 

iW^ = -U iF (A-27) 

c e 

The resulting Euler- Lagrange equation is 

AA- 

LF + CF = U (A-28) 

e 


Noting that the inertance and compliance are inverses of mass and stiffness, respec- 
tively, the natural frequency is 


U) = 



(A- 29) 


Example 2; 

Consider now the free-free system illustrated in Fig. A- 2. The kinetic and 
strain energies are 

1 Ao 1 A A 2 

U =-^Cf 2 (A- 30b) 

C 2 

Noting that the external impulse F^ is prescribed and has no variation (iFe = 0) the 
Euler- Lagrange equation is 

(L^ + Lg) F + CF = -Lg Fg (A-31) 


A- 7 



L=-“ « IIMERTANCE 



FREE BODY - 



KINETIC ENERGY: T 

c 2 


A 

F 


A 

F 


STRAIN ENERGY: U =4cF^ 

c 2 


1 


A - 

F, U VIRTUALWORK:6W^ = -U 5 
e C e 


Fig. A'1 Simple Oscillator with Base Excitation 
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Fig. A-2 Free-Free System 
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In terms of the particle masses, (Mj^, Mg)* the inertance of the system is 


(^1 



(A- 32) 


which is immediately recognized as the inverse of the relative mass. Rigid body 
motion singularities are automatically removed in the complementary formulation. 
The complementary forcing function 


A 1 A - 

L„F F =U 

2 e Mg e e 


(A-33) 


is equivalent to an effective base motion input. 

Example 3; 

Consider the two-degree of freedom system illustrated in Fig. A- 3. The kinetic 
and strain energies are 

(A -34a) 
(A -34b) 

(A- 35) 

The flexibility and inertance matrices are inverses of the stiffness and mass matrices 
respectively of the system when the relative displacements Xj and - Xj^) are taken 
as the generalized displacements. If in the complementary formulation the net 
impulses are taken as generalized coordinates, i.e., 

A A 

Fj"=F^ (A- 36a) 

A , A A 

^ 2 = Fg - F^ (A- 36b) 


1 ^91 ^ ^9 

^c = 2Vl "2 4(^2-^1> 


1 A O 1 A o 


and the Euler- Lagrange equations are 
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Fig. A-3 Supported Two-Degree of Freedom System 
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the kinetic and strain energies are expressed as 

1 2 1 ^2 

1 ^21 ^ A 2 

and the Euler- Lagrange equations are 



The above flexibility and inertance matrices are the inverses of the stiffness 
and mass matrices, respectively, of the system when the absolute displacements are 
taken as generalized coordinates. 

It is interesting to note that the use of absolute impulse as a generalized coor- 
dinate results in a relative displacement equivalent formulation while use of relative 
(net) impulse generalized coordinates results in an absolute displacement equivalent 
formulation. 
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APPENDIX B - HEAT CONDUCTION - FLUID FIX3W ANALOGY 

The current NASTRAN hydroelastic modeling capability is limited to axi- 
symmetrically shaped fluid configurations modeled with the cylindrical core and 
toroidal elements defined on CFLUID2, CFLUID3 and CFLUID4 connect cards. The 
possible need for a tilted liquid free surface fluid model for space shuttle external 
tank analysis necessitates the use of generally configured fluid finite elements. 
Tetrahedral and hexahedral finite elements would satisfy this requirement as well as 
any future requirement to model general non-axisymmetric fluid geometries. The 
desired general fluid elements are available in NASTRAN in the form of heat conduc- 
tion elements; these may be utilized as incompressible fluid dynamic finite elements 
by implementation of the analogy presented below. 

The field equation and flow boundary condition for an Incompressible, inviscid 
fluid undergoing small motion are; 

7 = 0 in V (B-1) 

and 

17 P . n = — = - pu on S (B-2) 

5 x^ ^ n 

A 

On the boundary surface (outward normal x^, n). The field equation and heat 
flux boundary condition for an isotropic solid under steady state thermal loading are: 


V^T= 0 

in V 

(B-3) 

dr 

i on S 

(B-4) 


n 


The analogy is now immediately apparent with the following variables taking 
equivalent roles: 

pressure, P: temperature, T 

outward normal acceleration, outward normal flux/area, f^ inverse 
mass density, p : inverse thermal conductivity, 1 /K.j. 
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Consider the "generalized energy" principles utilized to derive the finite 
element fluid flow and heat conduction equations. The fluid "generalized potential 
energy" is (Ref. 12): 

Uf =Y / ^(VP.rP)dV (B-5) 

V V 

the corresponding thermal potential function is (Ref. 12): 

Ut =J / K^(pr.7T)dV (B-6) 

■^V V 


in the special case of an isotropic material. 

The variation of outward flow potential for the fluid and the corresponding 
quantity for the thermal surface potential are; 

(5Uj =- /(JP(ipp).'SdS= f 6 


dS 


(B-7) 


B s 


and 


au =- /aT(KV'T)-ndS = /arf dS 


(B-8) 


B s 


respectively. 

The generalized energy expressions Eqs. B-5, B-6, B-7 and B-8 form the 
basis for derivation of mathematically equivalent incompressible fluid flow and 
thermal finite elements. The sets of algebraic equations describing fluid dynamic 
and thermal states, respectively, are 


K. P, = -S.U. 
f.. j i 1 
ij ^ 


and 


(B-9) 


T. = -S^f^ (B-10) 

ij •' 

The utility of fluid flow and thermal elements has just been proven to be inter- 
changeable; NARTRAN is capable of modeling general incompressible fluid flow con- 
figurations. The thermal elements required are defined by CTETRA, CHEXAl and 
CHEXA2 connect cards. Thermal material properties following the analogy (mass 
density, p : inverse of thermal conductivity, are specified on a MAT4 card. 
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